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The  development  of  a  wavelet-based  feature  extraction  technique  specifically  targeting  FOD- 
event  induced  vibration  signal  changes  in  gas  turbine  engines  is  described.  The  technique  performs 
wavelet  analysis  of  accelerometer  signals  from  specified  locations  on  the  engine  and  is  shown  to  be 
robust  in  the  presence  of  significant  process  and  sensor  noise.  It  is  envisioned  that  the  technique 
will  be  combined  with  Kalman  filter  thermal/health  parameter  estimation  for  FOD-event  detection 
via  information  fusion  from  these  (and  perhaps  other)  sources.  Due  to  the  lack  of  high-frequency 
FOD-event  test  data  in  the  open  literature,  a  reduced-order  turbofan  structural  model  (ROM)  was 
synthesized  from  a  finite  element  model  modal  analysis  to  support  the  investigation.  In  addition  to 
providing  test  data  for  algorithm  development,  the  ROM  is  used  to  determine  the  optimal  sensor 
location  for  FOD-event  detection.  In  the  presence  of  significant  noise,  precise  location  of  the  FOD 
event  in  time  was  obtained  using  the  developed  wavelet-based  feature. 

1.0  Introduction 


Ingestion  of  foreign  objects  by  a  turbofan  engine,  particularly  birds  and  ice,  is  one  of  the  leading 
causes  of  uncontained  rotor  events  in  commercial  jet  transport.  There  are  numerous  cases  of  bird 
ingestion  contributing  to  accidents,  some  of  them  fatal  [1].  Specific  procedures  have  been 
developed  for  situations  where  ingestion  is  suspected,  and  forensic  analysis  has  repeatedly  shown 
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that  they  would  have  been  the  appropriate  actions  to  take  in  cases  where  rejected  takeoffs 
motivated  by  bird  ingestion  resulted  in  accidents.  These  are  the  worst  situations;  in  the  vast 
majority  of  cases,  bird  ingestion  does  not  affect  the  safe  outcome  of  a  flight  and  may,  in  fact,  go 
unnoticed  by  the  flight  crew.  However  ingestion  does  pose  a  risk  of  Foreign  Object  Damage 
(FOD),  and  thus  is  important  to  identify,  if  possible.  When  a  multi-engine  aircraft  flies  through  a 
flock  of  birds,  potentially  damaging  more  than  one  engine,  the  pilot  needs  to  understand  the  status 
since  his  resulting  actions  may  be  different  depending  on  the  number  of  engines  involved.  Even  in 
a  case  where  there  is  no  apparent  damage  to  an  engine  as  observed  from  the  cockpit,  latent  effects 
(e.g.  cracks  that  can  be  propagated  by  high  cycle  fatigue)  may  be  present  that  require  maintenance 
action  [2]. 

The  critical  consequence  of  foreign  object  ingestion  is  engine  surge,  potentially  resulting  in  the 
loss  of  power.  The  flight  crew  can  recognize  foreign  object  ingestion  through  a  combination  of 
instrument  readings  and  sensory  cues.  These  include  such  symptoms  as  a  thud  or  bang,  a  fire 
warning,  a  visible  flame  coming  out  of  the  engine,  vibration,  yaw  of  the  airplane  caused  by  thrust 
imbalance,  high  Exhaust  Gas  Temperature  (EGT),  change  in  the  spool  speeds,  smoke/odor  in  cabin 
bleed  air,  and  Engine  Pressure  Ratio  (EPR)  change.  It  is  important  to  note  that  for  impact-type 
FOD  (due  to  ice,  birds,  runway  debris,  etc.),  the  damage  is  primarily  to  the  fan  and  front  part  of  the 
engine,  with  the  extent  of  the  damage  determined  by  the  geometry,  angle  of  impact,  hardness, 
relative  speed,  etc.  of  the  object. 

One  key  reason  for  detecting  FOD  automatically  is  that  most  such  events  occur  close  to  the 
ground  when  the  flight  crew  is  occupied  flying  the  plane.  Only  after  the  crew  stabilizes  the  aircraft 
at  a  safe  altitude  should  they  take  action  on  the  engine.  Once  the  aircraft  is  stable,  the  reduced 
workload  in  the  cockpit  environment  provides  better  circumstances  under  which  the  crew  can 
analyze  the  situation,  and  they  would  benefit  greatly  by  having  full  knowledge  of  the  engines 
involved  and  the  likelihood  of  the  event  based  on  data. 

Given  that  many  of  the  potential  symptoms  of  FOD  are  not  unique  to  that  type  event,  an 
automatic  system  for  FOD  detection  should  take  information  from  multiple  sources  to  provide 
confidence  in  its  diagnosis.  Just  like  a  pilot  does,  this  system  must  fuse  information  in  a  way  that 
provides  a  measure  of  the  likelihood  of  foreign  object  ingestion  in  order  to  determine  any 
corrective  action.  Additionally,  since  several  of  the  potential  symptoms  are  described  in  terms  of 
human  senses,  alternate  information  sources  need  to  be  developed.  These  sources  should  utilize  the 
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standard  engine  sensor  suite  (or  a  very  similar  suite)  to  address  issues  such  as  certification  and 
retrofit.  Naturally  processing  requirements  may  vary  from  current  systems  simply  due  to  the  fact 
that  additional  on-line  signal  processing  would  be  required,  but  the  impact  would  be  minimal  since 
such  a  system  would  not  be  flight  critical. 

For  approximately  two  decades  techniques  based  on  the  Kalman  filter  algorithm  have  been 
applied  to  turbofan  engine  diagnostics  [3,4].  Specifically,  analytical  (or  virtual)  measurements 
provided  by  Kalman  filter  estimators  have  been  used  for  detection  of  engine  degradation  via 
estimation  of  a  set  of  health  parameters,  which  are  generally  not  measurable  themselves  (e.g., 
compressor  efficiency)  and  are  calculated  via  knowledge  of  measurable  quantities.  The  degradation 
monitored  may  be  gradual  in  nature,  e.g.,  worn  components  that  increase  internal  clearances, 
resulting  in  decreased  component  efficiency.  It  may  also  be  abrupt  as  is  the  case  when  foreign 
objects  are  ingested  by  the  engine.  Changes  in  component  efficiencies,  high  and  low  spool  speeds, 
and  thermodynamic  parameters  have  been  determined/observed  during  foreign  object  damage 
(FOD)  events  [3,4].  In  some  cases,  these  parameter  changes  alone  may  not  be  conclusive  proof  that 
a  FOD  event  has  occurred.  Absent  from  these  investigations  is  the  incorporation  of  structural 
vibration  signals  as  a  means  to  aid  in  positive  identification  of  a  FOD  event. 

To  successfully  apply  available  sensor  fusion  techniques,  the  same  event  should  be  detected 
using  sensors  that  have  significantly  different  physical  characteristics  (e.g.,  thermocouples 
compared  to  accelerometers)  and  rely  on  measuring  physically  different  parameters  (e.g., 
temperature  and  acceleration)  [5,6].  It  would  appear  that  fusing  health  parameter  estimates  with 
structural  response  information  acquired  during  a  FOD  event  could  provide  conclusive  evidence 
that  a  FOD  event  had  occurred.  Furthermore,  fusion  of  virtual  measurements  and  vibration  signal 
features  may  enable  discrimination  between  a  normal  aircraft  maneuver/operating  condition  and  a 
more  severe  FOD  event  even  though  both  events  appear  to  have  similar  looking  health  parameter 
data  profiles  (due  to  inherent  limitations  of  the  diagnostic  model  used).  Similarly,  vibration  signals 
alone  may  possess  FOD-like  characteristics  depending  on  the  maneuver  e.g.,  a  rapid  change  in 
thrust  demanded  from  the  engine  would  result  in  transmission  of  temporary  impulse-like  force 
imbalances  through  the  engine  structure.  The  combination  of  traditional  health  parameter 
estimation  and  of  state-of-the-art  signal  processing  techniques  applied  to  structural  vibration 
signals  (e.g,  wavelet  analysis  [7,8])  could  provide  the  “finger  print”  necessary  to  positively  identify 
a  FOD  event. 
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This  paper  describes  the  development  of  a  wavelet-based  feature  (WF)  extraction  technique 
specifically  targeting  FOD-event  induced  vibration  signal  changes  in  gas  turbine  engines.  It  is 
envisioned  that  the  wavelet-based  feature  will  be  fused  with  measurements  or  estimates  of  engine 
thermal  parameters,  themselves  having  FOD-specific  response  characteristics,  to  provide  a  high- 
confidence  FOD  event  diagnosis.  Successful  automated  identification  of  FOD  events  would  be 
especially  useful  in  unmanned  (autonomous)  vehicle  applications,  where  the  absence  of  the  pilot 
precludes  the  use  of  human  monitoring  of  engine  conditions. 

To  the  authors’  knowledge,  engine-specific  vibration  signatures  associated  with  actual  FOD 
events  e.g.,  rotor  frequencies  excited  due  to  foreign  object  impact  on  the  fan  disk,  is  not  available 
in  the  open  literature.  Potential  sources  of  high-frequency  structural  data,  specifically  high-fidelity 
finite  element  models  of  aircraft  engines,  do  not  easily  lend  themselves  to  sensor  placement  and 
diagnostic  techniques  based  on  well  established  linear  system  theory.  As  well,  these  models  are 
typically  computation-time  intensive  due  to  the  vast  amounts  of  time-dependant  spatial  data  they 
produce  and  do  not  lend  themselves  for  direct  application  to  the  initial  stages  of  control  and 
diagnostic  (CDS)  development  and  testing.  To  facilitate  WF  development,  a  reduced  order  model 
is  synthesized  to  provide  the  structural  response  signals  indicative  of  a  FOD  event.  The  low-order 
structural  model  is  also  used  to  determine  the  optimal  accelerometer  locations  for  feature 
extraction.  This  requires  a  linear  model  in  state-space  form  to  facilitate  the  calculation  of  an 
Observability  Grammian  matrix  [9]  that  would  determine  whether  a  mode  excited  by  a  FOD  event 
could  be  detected  by  a  particular  sensor  arrangement. 

2.0  Development  of  reduced-order  turbofan  rotor  structural  models 

Finite  element  structural  analysis  codes  provide  invaluable  information  concerning  structural 
deformation,  stress  intensities  and  modal  behavior.  However,  finite  element  codes  tend  to  be 
cumbersome  and  time  consuming  to  work  with  when  used  to  design  and  test  control  and  diagnostic 
system  (CDS)  algorithms,  especially  when  performing  multiple  design  iterations.  State-space 
models,  on  the  other  hand,  are  the  form  of  choice  for  many  signal  processing  and  Multiple  Input 
Multiple  Output  (MIMO)  controller  applications.  A  method  for  reduced-order  modeling  of 
turbofan  rotor  structures  is  presented  which  uses  output  from  a  finite-element  code  in  the  form  of 
natural  frequencies  and  eigenvectors  to  create  a  state-space  model  suitable  for  use  in  CDS 


NASA/TM — 2004-2 13118 


4 


development  and  testing.  The  resulting  state-space  model  can  easily  be  ported  to  control  system 
development  software  such  as  the  MATLAB/SIMULINK®  suite  of  software  tools  [10].  The 
following  presents  the  reduced-order  model  development  process  as  applied  to  the  modeling  of  a 
generic  turbofan  rotor,  ultimately  used  for  generating  a  simulated  FOD  event  structural  response. 

For  the  purpose  of  developing  reduced-order  structural  models  for  application  to  CDS 
development,  a  dynamic  model  based  on  the  dominant  modes  of  vibration  of  a  turbofan  rotor¬ 
bearing  system  is  considered  adequate.  The  modal  behavior  characterized  by  the  model  should  be 
representative  of  that  presented  in  the  literature  for  turbofan  engines  [1 1,12].  For  the  present  study 
a  specialized  rotor  dynamics  finite  element  code,  DyRoBes  [13],  is  used  to  simulate  the  structural 
response  resulting  from  unbalance  and  foreign  object  damage.  This  code  is  suitable  for  the  present 
study  since  it  has  the  capability  to  model  critical  features  of  an  engine  system  and  can  sufficiently 
simulate  engine  dynamic  response,  particularly  the  response  at  sensor  locations.  Being  a  special- 
purpose  Finite  Element  code  i.e.,  designed  specifically  for  rotor  dynamic  analysis  and  users  with  a 
basic  knowledge  of  FE  methods,  it  is  also  relatively  straightforward  to  use. 

The  structural  model  used  for  this  study  was  derived  from  a  generic  engine  model  created  to 
develop  an  industry  standard  for  engine-airframe  structural  analysis.  The  model  is  of  a  dual  spool 
rotor  mounted  on  flexible  supports  (Figure  2.1).  The  flexible  supports  are  characterized  to  simulate 
the  static  engine  structure  and  wing. 

The  first  spool  is  representative  of  the  fan  and  low  pressure  turbine  stages  while  the  second 
spool  represents  the  high  pressure  compressor  and  turbine.  An  interstage  bearing  connects  the  two 
spools  together.  The  low  pressure  spool  is  1 12  inches  long  and  consists  of  18  shaft  finite  elements. 
A  large  disk  is  located  at  the  left  end  to  represent  the  fan  and  several  smaller  disks  are  located  on 
the  right  end  for  the  low  pressure  turbine.  The  entire  shaft  weighs  2494  pounds.  The  high  pressure 
spool  is  80  inches  long  and  consists  of  16  shaft  finite  elements.  There  are  eight  disks  on  the  shaft’s 
left  end  representative  of  a  compressor.  One  disk  is  located  on  the  shaft’s  right  end  representative 
of  the  high  pressure  turbine.  The  high  pressure  shaft  weighs  1280  pounds.  Each  shaft  is  supported 
on  bearings  as  shown  in  the  figure.  The  bearings  in  turn  are  supported  on  flexible  supports 
representative  of  the  support  flexibility  provides  by  the  static  engine  structure.  The  bearing  and 
support  stiffnesses  are  assumed  to  be  1.0  x  106  lb/in  and  1.0  x  103  Lb/in  respectively. 

The  transient  response  from  rotor  imbalance  and  FOD  of  the  turbofan  DyRoBes  model  is 
shown  in  Figure  2.2.  A  typical  rotor  vibration  sensor  signal  measured  during  a  FOD  event  (with  a 
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sensor  located  at  one  or  more  of  the  bearing  locations  shown  in  Figure  2.1)  would  consist  of  the 
response  of  a  dominant  low-frequency  lateral  mode  excited  from  an  impulsive  moment  at  the  fan 
disk  /  rotor  shaft  interface  (about  the  X-axis  in  Figure  2.3)  due  to  the  initial  impact  of  the  foreign 
object  upon  the  fan  disk.  Higher  frequency  lateral  modes  would  also  be  excited  due  to  the  high 
frequency  content  of  the  foreign  object  impact  and  would  appear  in  the  sensor 


Figure  2.1.  DyRoBes  FEM  schematic  of  turbofan  rotor-bearing  system 


signal  immediately  after  the  event.  Possible  longer  lasting  effects  due  to  resulting  permanent  blade 
damage  i.e.,  an  imbalance  force  occurring  at  a  frequency  corresponding  to  once  per  revolution  of 
the  rotor,  may  also  appear.  Rotating  equipment  tends  to  have  a  relatively  small  degree  of  inherent 
imbalance  after  manufacturing  [11],  which  will  also  be  included  in  the  reduced  order  structural 
model  for  a  “steady  state”  response.  To  the  authors’  knowledge,  high-frequency  accelerometer 
data  from  actual  turbofan  engine  FOD  events  does  not  appear  in  the  open  literature,  thus  the 
hypothesized  sequence  of  events  are  predominantly  based  on  the  authors’  experience  and  intuition. 
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The  following  section  describes  the  method  used  to  develop  a  Reduced  Order  Model  (ROM)  that 
closely  approximates  the  actual  rotor  response  but  is  more  amenable  to  real-time  assessment  of 
feature  extraction  techniques  used  for  identifying  FOD. 


Figure  2.2.  Example  FOD  event  FEM  model  vibration  response  characteristics.  Calculated 
displacements  are  functions  of  the  imbalance  magnitude  and  foreign  object  impact 
characteristics  (e.g.,  relative  speed,  size,  etc...). 

2.1  Dynamic  Representation  of  Rotor  Bearing  Systems  in  the  State  Space 


Dynamic  systems  may  be  represented  as  transfer  functions  or  in  state-space  [14].  For  CDS 
development  the  state-space  form  is  often  preferred  and  will  be  the  form  used  for  the  reduced  order 
model  synthesized  from  DyRoBes  FEM  modal  analysis  data.  For  a  multiple  degree  of  freedom 
(MDOF)  translational  system  the  equation  of  motion  in  physical  coordinates  is: 


[M]{q}+[c]{q}+[Klq}={F}  2.1 
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with  the  damping  matrix  C  typically  assumed  to  be  proportional  to  the  mass  and  stiffness  matrices 

[c]  =  a[M]+p[K]  2.2 

The  dynamics  of  rotating  systems  are  complicated  by  the  presence  of  gyroscopic  effects 
[11,15].  Expanding  equation  2.1  in  terms  of  translational  and  rotational  degrees  of  freedom  and 
including  gyroscopics  yields: 


M 

0 

0 

o' 

X 

^XX 

Cxy 

^xa 

CX0 

X 

0 

M 

0 

0 

y 

>  + 

Cyx 

cyy 

Cya 

Cy0 

y 

0 

0 

la 

0 

a 

Cax 

C<xy 

Caa 

ca0  “2IaQ 

a 

0 

0 

0 

ie_ 

0 

c0x 

C@y 

c0a  +2Iefi 

C00 

0 

* 

X 

X 

KXy 

KXa 

X 

CD 

X 

Fx 

Kyx 

Kyy 

Kya 

Kye 

y 

Fy . 

Kax 

Kay 

Ka« 

Kae 

a 

F« 

K0X 

o 

K0a 

K0e 

0 

Fe. 

Unlike  a  non-rotating  system,  the  damping  coefficient  matrix  in  Equation  2.3  is  a  function  of 
the  shaft  rotational  speed,  £2  i.e.,  gyroscopic  effects  result  in  the  rotor’s  lateral  natural  frequencies 
and  mode  shapes  being  a  function  Q .  To  determine  the  natural  frequencies  and  mode  shapes, 
assume  that  the  free  response  of  the  system  is  sinusoidal  in  nature  i.e., 

-  ,  r  Jo)  t 

{q(t)}={cp}e  n  2.4 

where  {<p}is  the  mode  shape  (eigenvector)  corresponding  tocon ,  the  natural  frequency  (eigenvalue). 
N  is  the  number  of  modes  included  in  the  model.  Inserting  Equation  2.4  into  the  unforced, 
undamped  version  of  Equation  2.3  gives 

(-(on2[M]+to„[G]+[K])[(p}=0  2.5 


NAS  A/TM— 2004-2 13118 


8 


where  the  gyroscopic  matrix  G  is  given  by 


0  0  0  0 

0  0  0  0 

0  0  0  -2ian 

o  o  2i en  o 


2.6 


A  trivial  solution  to  Equation  2.5  is  obtained  with{(p}=0.  For  a  nontrivial  solution 

(-<on2[M]+;io,,[G]+[K])=0  2.7 

must  be  satisfied.  To  obtain  a  reduced  order  linear  model  with  coefficients  that  are  independent 
of  Q. ,  i.e.,  with  constant  natural  frequencies coni  (i=l..N)  and  modal  damping,  it  is  assumed  that 
the  natural  frequencies  and  modes  shapes  calculated  via  Equation  2.7  at  the  rotor  critical  speeds 
adequately  account  for  the  effects  of  gyroscopics  over  the  operating  range  of  interest.  For  the  rotor 
described  in  Section  2.0,  the  coni  and  cp;  are  calculated  using  the  DyRoBes  FEM  code.  The  mode 
shapes  cpj  obtained  are  orthogonal  with  respect  to  the  mass  and  stiffness  matrices  and  provide  the 
following  normalization 


{<Pi  }T  [M]{(pi }  =  1  2.8 

foi}T[K]{(pi}=Wnj  2.9 

Given  the  assumptions  described  above,  the  equations  of  motion  for  the  linearized  rotational 
system  take  the  form  of  Equation  2.1.  For  lightly  damped  systems  the  matrix  C  is  typically 
assumed  to  be  diagonal  [16].  With  the  eigenvectors  mass-normalized  and  [<3>]  =  [<pj  •  •  *<pN  ], 
the  coordinate  transformation 

{q(t)M<b]{il(t)}  2.10 
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is  inserted  into  Equation  2.1  and  pre-multiplied  by  Ot  .  The  dynamics  of  the  system  in  modal 
coordinates  are 


{rj}+  diag[24jCOni  ]{fi}+  ding^  ]{q}=  Ot{f} 
with  the  corresponding  state-space  form 


n 


-diag 


2^i“ni 


11 

A. 


+<h 


iq}=[®  O^j 


and  the  modal  damping  ratios  4i  defined  as  [16] 


$■  = 


Cj 

C 

^  i  critical 


2.11 


2.12 


2.13 


For  the  present  study,  the  coordinates  of  interest  (Figure  2.3)  are 

*>>“{*} 

at  specified  points  (i.e.,  node  locations  in  Figure  2.1)  in  the  rotor-bearing  system. 

The  above  presents  a  simplified  formulation  for  the  equations  of  motion  of  a  rotational  system 
which,  in  a  linearized  sense,  include  the  gyroscopic  effects.  These  effects  are  accounted  for  in 
DyRoBes’  analyses.  Modal  damping  ratio  estimates  are  based  on  whirl  analysis  log-decrement 
calculations  at  the  rotor-bearing  system  critical  speeds  and  the  nominal  speed  of  rotation. 
Comparison  of  the  reduced-order  structural  model  frequency  response  (to  an  assumed  imbalance) 
and  transient  response  to  the  FEM-generated  data  will  provide  a  measure  of  the  fidelity  of  the 
reduced-order  model  to  appropriately  describe  a  FOD-event  induced  structural  response. 
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Figure  2.3.  Coordinate  system  used  in  the  FEM  and  reduced  order  state  space  structural  models. 
Cyclic  imbalance  forces  are  perpendicular  to  the  Z-axis.  Impulse  moments  due  to  a  FOD  event  are 
about  the  X-axis. 

2.2  Foreign  Object  Damage  Events  modeled  using  Reduced  Order  Rotor  Models 

from  FEM  output  data. 

The  state  space  model  presented  in  Equation  2.12  requires  the  modal  or  natural  frequencies  of 
the  system  and  the  eigenvectors  (i.e.,  modal  displacements)  at  specified  locations  on  the  rotor  or 
bearings  /  supports.  For  the  purposes  of  FOD  detection  algorithm  development,  the  locations  of 
interest  would  be  those  corresponding  to  possible  vibration  sensor  locations  i.e.,  nodes 
corresponding  to  the  bearings  and  supports  as  shown  in  Figure  2.1.  Commercial  FEM  codes  can 
readily  provide  this  information  [13].  Transformation  from  modal  coordinates  to  physical 
coordinates  is  performed  via  Equation  2.10,  with  the  matrix®  formed  by  columns  corresponding 
to  the  eigenvectors  for  each  frequency  being  considered.  To  transform  external  forces  on  selected 
locations  on  the  shaft  (or  supports)  from  physical  to  modal  coordinates,  the  external  force  vector  is 
pre-multiplied  by®T . 

Physical  (continuous)  systems  are  inherently  infinite  order.  Thus  FEM,  being  lower  order  than 
the  physical  system  being  modeled,  are  implicitly  inaccurate  with  FEM  being  more  accurate  in  the 
low  frequency  range.  The  frequencies  of  interest  are  from  0  to  20,000  rpm,  with  the  nominal 
rotational  speed  of  the  low  pressure  shaft  being  9000  rpm.  This  frequency  range  will  provide  the 
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required  response  due  to  fan  imbalance,  as  well  as  the  impulse  response  characteristics  due  to  FOD 
impact  upon  the  fan  disk.  For  adequate  imbalance  response  over  the  full  range  of  operating  speeds, 
the  natural  frequencies  and  mode  shapes  obtained  from  the  FEM  correspond  to  the  so-called  critical 
speeds  of  the  shaft  [  1 1],  To  attain  proper  response  to  a  FOD  event,  several  of  the  lateral  modes 
calculated  via  an  FEM  whirl  analysis  at  a  specified  rotor  speed  are  included  in  the  model  as  well. 
The  number  of  lateral  modes  included  is  determined  by  the  performance  of  the  ROM  with  respect 
to  the  FEM.  Once  again,  the  goal  of  the  present  study  is  to  develop  structural  models  to  test  real¬ 
time  diagnostic  algorithms,  with  these  models  providing  acceptable  response  characteristics  for  the 
forcing  functions  of  interest  (imbalance  and  foreign  object  impact)  while  maintaining  the  model 
order  low.  The  reduced-order  model  development  technique  described  above  can  be  applied  to 
FEM  developed  using  codes  intended  for  more  elaborate  analyses. 

The  reduced  order  dynamic  model  was  implemented  in  MATLAB®,  using  the  companion 
software  package  SIMULINK®,  in  state-space  form.  Figure  2.4  depicts  the  SIMULINK® 
structural  model.  Gas  turbine  rotors  tend  to  have  an  inherently  small  degree  of  imbalance,  which  is 
modeled  as  an  input  corresponding  to  a  variable  amplitude  sine  wave  with  a  frequency  equal  to  the 
high  and  /  or  low  speed  rotor  speeds.  Imbalance  force  amplitude  is  frequency  dependent  according 
to  the  following  relation 


^imbalance  =  ^  Q 


2.15 


where  m  is  the  total  mass  of  the  rotating  component,  e  is  the  effective  eccentricity  in  inches,  and  Q. 
is  the  shaft  rotational  speed  in  rad-s-1.  If  a  FOD  event  results  in  additional  imbalance,  the 
eccentricity  is  increased  by  an  additional  amount.  For  the  present  study,  the  fan  is  assumed  to  have 
an  eccentricity  (i.e.,  inherent  imbalance)  of  0.001  inch. 

A  FOD  event  is  modeled  as  an  object  hitting  a  point  along  the  radius  of  the  fan,  resulting  in  a 
moment  input  at  the  node  on  the  rotor  corresponding  to  the  fan  location  (Node  1  in  Figure  2.1).  The 
magnitude  of  the  input  moment  is  determined  as  follows.  Given  a  foreign  object  of  mass  oifo  with 
an  initial  linear  velocity  relative  to  the  fan  disk  vFo  along  the  z-axis  (Figure  2.3),  conservation  of 
angular  momentum  may  be  applied  to  determine  an  estimate  of  the  initial  angular  velocity©  of  the 
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Figure  2.4.  MATLAB/SIMULINK  implementation  of  the  reduced  order  turbofan 
rotor/bearing  system  structural  model. 
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fan  disk  about  an  axis  perpendicular  to  the  axis  of  rotation  of  the  rotor  (X-axis  in  Figure  2.3)  due  to 
the  impact  i.e., 

e  =  rm  111 -ill  2.16 


where  I0  is  the  moment  of  inertia  of  the  fan  disk  about  the  X  axis  (Figure  2.3)  and  r  is  the  distance 
of  the  impact  from  the  axis  of  rotation  (Z  axis  in  Figure  2.3).  The  corresponding  moment  applied  to 
the  rotor  at  node  1  from  the  impulse  load  on  the  fan  disk  is  determined  by  the  rate  of  change  of 
angular  momentum,  given  a  specified  time  interval  for  the  impact  to  occur  At . 


F*=I* 


A0 

At 


2.17 


The  above  treatment  provides  a  straight-forward  means  for  determining  the  impulse-load  related 
torque  input  to  the  rotor  during  the  initial  stages  of  the  FOD  event,  and  is  considered  to  be  adequate 
within  the  level  of  complexity  of  the  model  needed  for  real-time  diagnostic  system  development. 

As  mentioned  previously,  information  regarding  specific  FOD  events  does  not  appear  in  the  open 
literature.  The  FOD  event  impact  time  interval,  At ,  is  dependant  upon  random  factors  such  as 
foreign  object  geometry,  location  and  angle  of  impact,  foreign  object  hardness,  relative  speed,  etc. 
and  is  estimated  based  on  the  authors’  “intuitive  feel”  for  an  event  that  would  cause  damage. 

Figure  2.5  presents  a  frequency  response  comparison  of  4,  10,  and  25  degree  of  freedom 
(DOF)  reduced  order  models  with  the  frequency  response  of  the  154  DOF  DyRoBes  FEM  (refer 
to  Figure  2.1  for  node  locations).  Figures  2.6  and  2.7  present  the  displacement  response  of  the 
10  DOF  ROM  and  FEM  to  a  combined  FOD  event  -  imbalance  input.  Figure  2.6  presents  the 
response  at  bearing  /  support  location  37  to  a  “large”  FOD  event  (mFo  =  2  lbm)  and  a  “small”  FOD 
event  (mFO  =  0.25  lbm).  Figure  2.7  presents  the  response  at  bearing  /  support  location  39  to  a 
“large”  FOD  event  (mFO  =  2  lbm)  and  a  “small”  FOD  event  (mFO  =  0.25  lbm).  The  imbalance  input 
is  calculated  using  Equation  2.15  with  an  eccentricity  of  0.001  inches  concentrated  at  the  fan  node 
(node  1  in  Figure  2.1).  The  magnitude  of  the  FOD  event  moment  input  at  node  1  is  calculated  using 
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Equation  2.16  with  a  radius  of  20  inches,  specified  mass,  a  At  of  0.1  seconds,  and  a  relative 
velocity  of  300  mph.  To  provide  adequate  FOD  event  response  characteristics,  the  first  lateral 
mode  from  an  FEM  whirl  analysis  at  9000  rpm  was  included.  The  10  DOF  model  provides 
acceptable  correspondence  to  the  FEM  and  is  considered  sufficient  for  the  present  study. 
Differences  between  the  reduced  order  model  and  the  FEM  are  due  to  the  assumptions  made,  i.e., 


Figure  2.5.  Steady  state  imbalance  response  comparison  of  4,  10,  and  25  degree 
of  freedom  (DOF)  reduced  order  models  for  two  different  bearing/support 
locations.  Imbalance  has  an  eccentricity  of  0.001  inch  located  at  the  fan. 
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(mFO  =  2  lbm),  and  (b)  a  “small”  FOD  event  (mFo  =  0.25  lbm). 
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Figure  2.7.  Response  at  bearing  /  support  location  39  to  (a)  a  “large”  FOD  event 
(mFo  =  2  Ibm)  and  (b)  a  “small”  FOD  event  (mFo  =  0-25  lbm)  (b). 
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the  ROM  uses  constant  modal  damping  ratios  and  natural  frequencies  over  the  full  range  of 
operating  speeds. 

For  both  locations,  the  change  in  the  vibration  signal  due  to  the  large  FOD  event  is  easily 
discemable  from  the  nominal  vibration  signal  due  to  inherent  imbalance  alone.  For  the  small  FOD 
events,  there  is  less  of  a  difference  from  the  normal  operating  characteristic.  For  smaller  FOD 
events  with  the  presence  of  sensor  or  process  noise  (which  would  be  the  case  in  an  actual  engine) 
the  event  could  easily  be  masked.  Thus,  the  signal  analysis  technique  employed  to  detect  such  an 
event  must  be  capable  of  identifying  the  characteristic  vibration  signature  in  the  presence  of 
process  and  sensor  noise. 

3.0  Optimal  Sensor  Placement  for  FOD  event  detection. 

Signal  processing  techniques  that  provide  information  concerning  the  internal  modes  (i.e., 
states  )  of  a  system  require  that  these  state  variables  be  observable  using  the  available  sensor  suite 
[17].  For  FOD-event  detection,  the  vibration  sensor  (e.g.,  accelerometer)  signal  must  contain  the 
modal  response  of  the  rotor  in  the  frequency  range  of  interest.  The  modes  of  a  system  are 
observable  if  and  only  if  the  matrix 

C 

CA 
CA"-1 

known  as  the  observabilty  matrix,  has  rank  n,  where  n  is  the  dimension  of  the  state  space  of  the 
linear  time  invariant  system  under  consideration.  For  the  10  mode  ROM  representing  the  rotor¬ 
bearing  system  shown  in  Figure  2.1  (with  nodes  37  through  41  as  candidate  accelerometer 
locations)  the  system  is  completely  observable.  To  determine  the  optimal  (or  worst  case) 
accelerometer  location  for  FOD  event  detection,  a  measure  of  the  degree  of  observability  of  a 
system  for  a  given  sensor  arrangement  must  be  determined.  One  method  suggested  in  [17]  requires 
calculation  of  the  Observability  Grammian  to  determine  which  states,  starting  from  a  specified  set 
of  initial  conditions,  have  an  influence  on  the  output  energy  E(x(0))  of  the  system  over  all  time. 
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Consider  the  following  measure  of  the  output  energy  at  some  initial  condition  x(0)  for  the  system 
of  Equation  2.3  under  observation  with  no  external  input 

00^00 

E(x(0))  =  jj|y(t)||  dt  =  Jy(t)Ty(t)dt  =  x(0)TQ  x(0)  3.2 

0  0 

where  the  Observability  Grammian  is  defined  as 

OO 

Q  =  JeAT‘CTCeAt  dt  3.3 

0 

For  a  stable  system,  the  existence  of  Equation  3.3  is  guaranteed  and  Q  may  be  found  by  solving 
the  following  Lyapunov  equation 


AtQ  +  QA  +  CtC  =  0  3.4 

For  a  reduced  number  of  accelerometers  (reduced  relative  to  the  initial  set  i.e.,  accelerometers 
37  through  41  in  Figure  2.1)  the  system  is  observable  if  and  only  if  Q  is  positive  definite.  With 
observability  established  for  the  reduced  sensor  suite  system,  Equation  3.2  is  used  to  determine  the 
output  energy  contribution  of  each  state  variable.  Large  output  energies  correspond  to  easily 
observable  state  variables  (modes).  The  optimal  suite  of  sensors  would  naturally  be  that  with  the 
(relatively)  largest  output  energies  for  the  modes  needed  to  identify  a  FOD  event.  Consider  the 
calculation  of  Equation  3.2  for  locations  37  and  40.  With  the  entry  in  the  initial  state  vector  (i.e., 
x(0))  for  a  specified  mode  set  to  1 .0  and  the  initial  conditions  for  each  of  the  other  modes  set  to  0.0, 
E(x(0))  for  each  mode  observed  separately  is  simply  the  diagonal  entry  of  Q  corresponding  to  each 
mode.  Thus,  a  judgment  as  to  the  optimal  placement  of  the  accelerometers  may  be  made  based  on 
the  diagonal  entries  of  Q  alone.  For  locations  37  and  40 
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diag[Q] 
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4.121  le- 003 
5.1742e-  003 
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4.3813e-004 
1 .6479e  -  004 
2.9526e  -  004 
2.2466e  -  008 
1.0822e-004 


diag[Q] 

Node  40 


3.5 


which,  due  to  the  higher  observation  energies  associated  with  the  first  three  modes,  implies  that 
placement  of  an  accelerometer  at  location  37  may  be  more  appropriate  for  detecting  a  FOD  event. 
Upon  performing  the  calculation  of  Equation  3.2  for  nodes  37  through  41  as  candidate 
accelerometer  locations,  node  37  was  determined  to  be  the  preferred  location  with  node  40  being 
the  least  preferred.  Thus,  the  analysis  (i.e.,  feature  extraction)  technique  selected  should  be  applied 
to  the  signal  generated  by  “accelerometer  37”  in  order  to  have  the  highest  likelihood  of  detecting  a 
FOD  event.  However,  placement  of  a  sensor  at  the  optimal  location  may  not  be  physically 
attainable,  thus  it  is  desirable  that  this  same  technique  should  be  robust  enough  to  detect  the  event 
at  a  suboptimal  (possibly  the  most  limiting)  location.  For  some  systems  this  conclusion  may  be 
drawn  directly  from  physical  intuition.  For  more  complex  systems  (which  may  not  easily  lend 
themselves  to  physical  intuition)  the  analytical  treatment  presented  above  provides  a  quantitative 
measure  for  optimal  sensor  placement. 


4.0  FOD  Event  Detection  Via  Rotor  Bearing  Vibration  Signal  Analysis 

Consider  the  simulated  bearing  accelerometer  signals  shown  in  Figure  4.1.  To  demonstrate  the 
robustness  of  the  WF  extraction  technique,  the  analyses  described  below  were  performed  on  the 
signal  from  the  most  limiting  accelerometer  location  on  the  rotor  bearing  system  considered,  which 
was  determined  from  the  sensor  placement  analysis  of  section  3.0  to  be  at  node  40  (Figure  2.1). 
This  location  would  result  in  the  weakest  signal  and  necessarily  correspond  to  the  most  difficult 


NAS  A/TM— 2004-2 13118 


20 


location  for  FOD  event  detection.  At  3.44  seconds,  a  0.5  lbm  foreign  object  hits  the  fan  disk  with  a 
speed  of  300  mph  (parallel  to  the  Z-axis  in  Figure  2.3)  at  a  radius  of  20  inches.  The  event  is 
modeled  as  a  pulse  of  width  0.04  seconds,  with  a  magnitude  determined  using  the  method 
presented  in  Section  2.2.  A  fan  disk  eccentricity  of  0.001  inches  is  assumed.  As  shown  in 
Figure  4.1a  ,  the  event  is  barely  noticeable  in  the  time  trace  for  a  noise-free  situation.  However 
accelerometer  signals  on  in-service  engines  are  typically  noisy.  In  addition  to  “process  noise”  i.e., 
random  motion  of  the  aircraft  due  to  wind  gusts  and  compensating  maneuvers  being  transmitted 
from  the  airframe  to  the  engine,  there  is  a  significant  amount  of  sensor  noise  which  may  mask  the 
occurrence  of  a  FOD  event.  Figure  4.1b  shows  the  accelerometer  signal  at  location  40,  with  a 
signal-to-noise  ratio  of  approximately  3.5  (1 1  dB).  Observation  of  the  signal  shows  no  well-defined 
point  in  time  at  which  one  would  identify  a  FOD  event  occurrence. 


01 _ i _ i - 1 - 1 - >— 

3.3  3.35  3.4  3.45  3.5  3.55 

Time  (sec) 


Figure  4.1 .  Bearing  40  accelerometer  output  signal  corresponding  to  a  FOD  event, 
a)  Noise  free  system  response,  b)  System  response  with  process  and  sensor  noise. 
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4.1  FOD  event  Detection  using  Fourier  Analysis  of  Rotor  Bearing  Vibration  Signals 


Vibration  signal  analysis  is  an  invaluable  tool  for  identifying  mechanical  component  problems. 
The  most  utilized  tool  in  vibration  analysis  i.e.,  the  Fourier  transform,  essentially  assumes  that  the 
signals  being  analyzed  are  of  infinite  duration  [18,19].  When  the  continuous-time  Fourier  transform 

Y(F)=  f  y(t)e~j2;rFtdt  4.1 


or,  more  practically,  its  discrete-time  counterpart 


Y(f)=  X  y(n)e“j2*fn 

n=-oo 


4.2 


is  applied  to  a  signal  of  finite  duration,  spectral  leakage  effects  can  significantly  reduce  the 
resolution  of  the  corresponding  power  (frequency)  spectrum  given  by 


E  y(n)e-J2*f" 


2 


n=  —  <*> 


4.3 


Windowing  of  the  original  signal  aids  in  mitigating  the  effects  of  finite  signal  length  (spectral 
leakage),  however  if  the  component  of  the  signal  being  analyzed  has  “compact  support”  (e.g.,  an 
impulse  due  to  a  FOD-event  induced  structural  response),  windowing  of  the  signal  may  not 
significantly  aid  in  identifying  the  occurrence  of  the  event  in  time  due  to  the  Fourier  transform 
basis  functions  i.e.,  sinusoids,  being  themselves  functions  of  infinite  support  and  limited  in  their 
ability  to  decompose  such  a  finite-duration  signal.  The  short-time  Fourier  transform  (and  the 
corresponding  time-dependant  power  spectrum)  uses  a  windowed  version  of  y(n)  to  compute  Y(f ) 
in  an  attempt  to  localize  a  change  in  a  signal  over  time  [19].  The  time-dependant  power  spectrum 
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for  the  noise-free  signal  (Figure  4.2  a)  shows  a  significant  change  at  approximately  3.3  seconds, 
which  disappears  as  time  progresses.  The  time-dependant  power  spectrum  of  the  noisy  signal 
(Figure  4.2  b)  shows  no  noticeable  change  in  time,  with  the  effect  of  the  FOD  event  “buried”  in  the 
spectra  due  to  the  noise.  Thus,  there  is  a  need  for  an  alternative  signal  analysis  technique  for  FOD 
event  detection. 


Peak  corresponding  to 


Figure  4.2.  Power  spectra  of  bearing  40  accelerometer  output  signal  corresponding  to  a 
FOD  event,  a)  Noise  free  system  response,  b)  System  response  with  process  and  sensor 
noise. 


NAS  A/TM— 2004-21 3 1 1 8 


23 


4.2  FOD  event  Detection  using  Wavelet  Analysis  of  Rotor  Bearing  Vibration  Signals 


Over  the  past  decade  the  discrete  time  wavelet  transform  (DTWT)  has  been  applied  to  a  wide 
range  of  signal  analysis  problems  e.g.,  de-noising  of  signals  as  well  as  time  localization  and 
reconstruction  of  short  duration  changes  [8].  Because  of  these  characteristics,  the  DTWT  is 
considered  to  be  a  viable  candidate  for  identifying  a  short  duration  change  in  a  bearing  vibration 
signal  due  to  a  FOD  event  corrupted  by  noise. 

The  DTWT  is  the  discrete-time  counterpart  of  the  continuous-time  wavelet  transform 

CWT(a,  r)  =  -j=  Jy(t)  4^-^jdt  4.4 

where  T'(t)  is  the  basic  or  mother  wavelet  and  ¥((t  -  x)/a)/Va  are  the  wavelet  basis  functions  which 
are  scaled  via  the  parameter  a  (by  compressing  or  stretching  the  mother  wavelet)  and  shifted  in 
time  by  x .  The  operation  performed  in  equation  4.4  on  the  original  signal  y(t)  may  be  interpreted 
as  the  cross  correlation  of  the  signal  y(t)  with  ^(t/a )/Va  ,  shifted  by  x/a  [18].  Thus,  the  CWT 

computes  the  component  of  y(t)  that  is  similar  to  vF(t/a)/Va  .  If  little  (or  no)  similarity  exists 
between  the  two,  then  the  CWT  will  be  small  (or  zero).  Larger  values  of  CWT  indicate  better 
correlation.  Figure  4.3  presents  an  overview  of  the  Wavelet  Transform  process  with  multiple  scale 
decomposition.  The  wavelet  family  used  is  the  Daubechies  family. 

Typically,  the  resulting  data  is  presented  in  terms  of  time  (or  k  if  in  discrete  time)  and  scale 
(i.e.,  coefficient)  of  the  chosen  wavelet  transform.  “Pseudo”  frequency  components  of  y(t), 
corresponding  to  the  scale  of  the  constituent  wavelets,  Fm  ,are  obtained  by  the  following 
expression  [20] 


TF; 

m 


4.5 
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where  m  is  the  scale,  T  is  the  sampling  interval,  and  Fc  is  the  wavelet  center  frequency  calculated 

by  taking  the  maximum  of  the  modulus  of  the  wavelet’s  Fourier  transform. 

Most  of  the  functions  for  'F(t)  have  been  developed  over  the  past  two  decades  with  the  Haar 

wavelet  being  the  first  documented  function  to  be  used  in  this  context[19].  Choice  of  'F(t)  is 

heavily  dependant  on  the  nature  of  the  signal  (or  a  localized  feature  of  a  signal)  being  analyzed.  For 

example,  if  a  vibration  signal  were  to  be  analyzed  to  determine  whether  a  FOD  event  had  occurred, 

one  would  expect  to  see  an  impulse-like  localized  feature  embedded  in  a  normally  noisy  sinusoidal 

signal.  Intuitively,  one  would  choose  a  wavelet  that  had  an  impulse-like  character  as  opposed  to 

one  that  resembled  a  step  function,  as  is  the  case  with  the  Haar  wavelet  mentioned  previously. 

Inverse  Wavelet  transform 

gives  signal  decomposition  in  terms 

of  time  and  scale  (pseudo  frequency) 


Figure  4.3  Example  bearing  vibration  signature  and  corresponding  wavelet 
decomposition  accompanying  a  FOD  event 

In  discrete-time  the  wavelet  transform  takes  the  form 

_m 

DTWTD(m,  n)  =  a0  2>(k)¥(a0-mk-nT0)  4.6 

k 
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where  m  is  the  scale  (degree  of  dilation)  and  n  corresponds  to  the  net  translation  in  time  (in  terms 
of  the  number  of  time  x0 )  of  the  wavelet  at  a  specific  scale.  For  scale  m=0,  the  coefficients  of  the 
so-called  scaling  function,  from  which  all  of  the  wavelets  in  a  given  family  are  related,  are 


DTWTA(0,  n)  =  £  y(k)cp(k  -  nx0 )  4.7 

k 

cp(k  -  nT0)is  typically  determined  via  iteration  using  the  following  dilation  equation 


N-l 

<p(k)  =  Xcn(P(2k  -  nx0) 


n=0 


4.8 


where  N  is  the  number  of  terms  in  cp(k)  .  The  scaling  function  cp(k)  determines  the  character  of 
the  wavelet.  The  characteristic  of  interest  in  the  signal  being  analyzed  typically  dictates  the  form  of 
the  wavelets  used.  A  wavelet  \|/(k)  is  related  to  it’s  associated  (p(k)  via 


\|/(k)  =  ^(-l)ncn(p(2k  +  nT0  -N  +  l)  4.9 

n=0 

Equation  4.9  is  also  a  dilation  equation  with \p(k) referred  to  as  a  dilation  wavelet  [19].  The 
scaling  function,  it’s  dilates,  the  corresponding  wavelets  and  their  dilates  are  orthogonal  to  one 
another  and  therefore  constitute  a  set  of  functions  by  which  an  arbitrary  function  may  be  built  [19], 
completely  analogous  to  the  Fourier  transform  and  series.  The  base  value  of  scaling  parameter,  a0 , 

and  the  time  shift,  x0 ,  are  typically  set  equal  to  2  and  1  respectively  for  computational  efficiency, 

analogous  to  the  Fast  Fourier  Transform  (FFT). 

The  convolution  of  Equations  4.6  and  4.7  is  commutative  i.e., 
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4.10 


DTWTA(0,n)  =  ^(p(k)  y(k  -  nx0) 
k 

m 

DTWTD(m,n)  =  a0“72>(k)y(a0-mk  -nx0)  4.11 

k 

Thus,  the  discrete-time  wavelet  transform  may  be  viewed  as  a  filtering  operation  with 
cp(k)  containing  the  coefficients  of  a  low-pass  Finite  Impulse  Response  (FIR)  filter  and  \|/(k)  the 
coefficients  of  a  high-pass  FIR  filter.  The  coefficients  corresponding  to  the  low-frequency  portion 
of  the  input  signal  i.e.,  approximation  coefficients  (Equation  4.10),  and  the  high-frequency  detail 
coefficients  (Equation  4.1 1)  at  a  given  scale  in  a  wavelet  decomposition  are  provided  by  successive 
low  pass  and  high  pass  filtering  operations  respectively.  Reconstruction  of  the  approximations  and 
details  at  a  specified  scale  m  is  performed  via  the  corresponding  inverse  wavelet  transform 


Dm(k)  =  XDTWTD(m,n)xivn(k) 

n 

4.12 

Am  (k)  =  £DTWTA(m,n)  ^  (k) 

4.13 

n 


and  is  implemented  via  banks  of  reconstruction  or  synthesis  filters.  Figure  4.4  presents  a  2-scale 
DTWT-based  signal  analysis  and  synthesis  using  an  FIR  quadrature  mirror  filter  bank.  Locations  in 
the  filter  bank  correspond  to  the  signal  characteristics  of  interest  (details  or  approximations). 
Scaling  is  accomplished  by  downsampling  after  high  or  low  pass  filtering.  Successive  upsampling 
and  filtering  performs  the  reconstruction.  When  the  reconstruction  filter  coefficients  correspond  to 
the  quadrature  mirror  filter  of  the  decomposition  filter,  perfect  signal  reconstruction  is  guaranteed 
[21].  This  is  the  practical  implementation  of  the  discrete-time  wavelet  transform  algorithm 
developed  by  Mallat  [22,  23,  24].  The  information  gained  from  the  decomposition  depends  on  the 
application  e.g.,  a  signal  may  be  “denoised”  by  omitting  the  detail  from  one  or  more  scales  during 
the  reconstruction.  As  well,  abrupt  changes  in  the  signal  e.g.,  a  FOD  event-induced  accelerometer 
response,  might  be  detected  by  observing  one  or  more  of  the  approximations. 
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Figure  4.5  presents  a  wavelet  analysis  of  the  signal  analyzed  in  Section  4.1,  the  noise-free 
bearing  accelerometer  signal  at  location  40.  The  Daubechies  8  wavelet  was  chosen  (a  subjective 
choice)  for  the  analysis  due  to  its  correspondence  to  the  signal  characteristic  of  interest.  Observing 
the  approximation  at  scale  8  (Figure  4.5b)  demonstrates  the  wavelet’s  ability  to  determine  the 
location  in  time  of  the  event.  Several  other  approximations  and  details  at  various  scales  (not  shown) 
also  revealed  the  occurrence  of  a  FOD  event.  The  analysis  of  the  noisy  signal  at  the  same  scale’s 
approximation  (Figure  4.6),  however,  illustrates  that  the  wavelet  chosen  (as  well  as  several  other 
wavelets  tested,  not  presented  here),  may  have  difficulty  highlighting  a  subtle  change  buried  in  a 
low  signal-to-low  noise  ratio  (SNR  =  3.5)  signal.  Additional  “conditioning”  of  the  decomposed 
signal  is  required. 


Scale  2  Approximation 


Low  Pass  Reconstruction  FIR 
Filter 


Low  Pass  Decomposition  FIR 
Filter 


~~\  High  pass  Reconstruction  FIR 
HK  Filter 


High  Pass  Decomposition  FIR 
Filter 


Downsample 
by  a  Factor  of  2 


Upsample 
by  a  Factor  of  2 


Figure  4.4.  Example  of  a  2-scale  DTWT-based  signal  analysis  and  synthesis  using  an  FIR 
quadrature  mirror  filter  bank. 
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Figure  4.5.  Wavelet  analysis  of  noise-free  bearing  accelerometer  signal  at  location  40. 
a)  accelerometer  signal,  b)  Corresponding  Daubechies  8  wavelet  inverse  transform,  scale  8 
approximation. 
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Figure  4.6.  Wavelet  analysis  of  noisy  bearing  accelerometer  signal  at  location  40. 

a)  accelerometer  signal,  b)  Corresponding  Daubechies  8  wavelet  inverse  transform,  scale  8 

approximation. 
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Abrupt  changes  in  an  otherwise  smooth  signal  will  correspond  to  a  localized  change  in  the 
signal  derivative  with  respect  to  time  (or  some  space-related  dependant  variable  depending  on  the 
application).  The  more  abrupt  the  change,  the  larger  the  local  value  of  the  derivative  with  respect  to 
time.  Additionally,  for  the  same  change  the  value  of  higher  order  derivatives 


dny(t) 

dtn 


4.14 


will  provide  a  frequency-dependant  gain  i.e.,  the  more  abrupt  the  signal  change,  the  larger  the  gain. 
This  quality  of  the  derivative  is  more  easily  observed  in  the  frequency  domain  by  taking  the 
Laplace  transform  of  Equation  4. 14  and  setting  s  =  jco 

=(jco)ny(jco)  415 

S  =  jft) 

The  corresponding  frequency-dependant  magnitude  is 


L[dny(t) 
[  dtn 


MAG^(jco)ny(jco)j  =  (On  MAG(y(  jco)) 


4.16 


Thus,  multiple  differentiation  in  the  time  domain  provides  an  “s11”  multiplication  factor  in  the 
frequency  domain.  The  higher  frequency  characteristics  of  a  signal  (e.g.,  a  FOD-event  would  be 
classified  as  a  high-frequency  event)  would  have  a  significantly  higher  gain  relative  to  the  low- 
frequency  components.  It  would  seem  that  in  a  noise-free  signal,  the  occurrence  of  a  FOD-event 
resulting  in  an  impulsive  input  to  the  turbofan  rotor  of  the  present  study  could  easily  be  identified 
using  this  technique.  Multiple  differentiation  of  noisy  signals  is  typically  avoided  in  most  signal 
processing  applications  due  to  the  amplification  provided  by  Equation  4.16.  Indeed,  multiple 
differentiation  of  the  signal  shown  in  Figure  4.6a  would  make  identification  of  a  FOD-event 
virtually  impossible,  primarily  due  to  amplification  of  the  noise. 
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Mallat  [25]  proposed  a  technique  for  identifying  a  wavelet  specifically  designed  for  edge 
detection  in  computer  vision  systems.  A  signal  (or  2-D  image  in  the  original  application)  is  passed 
through  a  smoothing  function  (e.g.,  a  low  pass  filter)  and  differentiated  multiple  times.  The 
combination  of  smoothing  and  differentiation  results  in  a  wavelet  customized  to  the  application. 
For  this  research,  the  technique  was  adapted  to  identify  the  FOD  event  induced  short-time  change 
in  accelerometer  signals.  The  smoothing  filter  chosen  for  the  present  investigation  is  the 
Daubechies-8  wavelet  decomposition.  As  shown  in  Figure  4.7,  the  accelerometer  signal  is  passed 
through  an  approximation  filter  bank  three  times  and  subsequently  passed  through  a  detail  filter. 
This  provides 


Low  Pass  Reconstruction 
FIR  Filter 


Low  Pass  Decomposition 
FIR  Filter 


jjD  I  High  Pass  Decomposition 
- '  FIR  Filter 


Downsample 
by  a  Factor  of  2 


Upsample 
a  Factor  of  2 


Figure  4.7.  Wavelet  transform-based  vibration  feature  extraction  implemented  using 
an  analysis  and  synthesis  FIR  filter  bank. 
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a  means  to  focus  on  a  subband  of  the  original  signal  where  the  feature  of  interest  is  thought  to  lie, 
the  net  effect  being  a  low-pass  filtering  or  smoothing  operation.  The  result  is  differentiated  three 
times  to  produce  the  wavelet-based  feature  of  interest.  Figure  4.8  presents  the  wavelet  feature 
extracted  from  the  accelerometer  signal  shown  in  Figure  4.6,  using  the  custom  wavelet  described 
above.  The  feature  extracted  is  the  spike  shown  in  Figure  4.6b.  The  exact  time  of  the  event  is 
precisely  determined  using  this  technique.  Indeed,  the  event  is  modeled  by  imposing  a  pulse,  whose 
magnitude  is  a  function  the  of  foreign  object  impact  characteristics  mentioned  previously,  on  the 
fan  disk.  The  two  spikes  shown  in  Figure  4.8  result  from  the  rising  and  descending  parts  of  that 
pulse  “filtered”  through  the  mechanical  (i.e.,  rotor-bearing)  system.  The  event  input  profile  to  the 
model  (i.e.,  the  pulse  input  described  above)  was  considered  to  be  appropriate  for  a  significant, 
structurally  damaging  event  i.e.,  where  the  foreign  object  absorbs  virtually  none  of  the  energy  of 
impact,  resulting  in  complete  and  immediate  transfer  of  energy  to  the  fan  disk.  This  would,  for 
example,  be  the  case  for  hard  objects  such  as  ice.  For  other  scenarios,  pulses  of  lower  magnitude  or 
with  non-infinite  rising  (or  descending)  slopes  may  require  higher-order  differentiation  for 
detection. 

In  summary,  the  result  of  the  following  sequence  of  operations: 

•  wavelet  decomposition  of  a  noisy  ( e.g.,  bearing  accelerometer)  signal 

•  reconstruction  of  selected  components  (approximations  or  details  at  a  selected  scales)  of  the 
signal 

•  multiple  differentiation  of  the  reconstructed  components 

can  precisely  time-locate  an  abrupt,  high-frequency,  change  in  the  signal.  However,  choice  of  too 
high  a  scale  may  result  in  the  high-frequency  event  being  filtered  out  of  the  signal. 
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Figure  4.8.  Wavelet  feature  extracted  from  noisy  bearing  accelerometer  signal  at  location 
40.  a)  accelerometer  signal,  b)  Corresponding  conditioned  Daubechies  8  wavelet  inverse 
transform,  scale  8  approximation. 
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5.0  Conclusion 

This  paper  describes  the  development  of  a  wavelet-based  feature  (WF)  specifically  designed  to 
identify  the  occurrence  of  FOD-events  in  gas  turbine  engines.  The  WF  is  extracted  from  rotor 
bearing  accelerometer  signals  and  is  shown  to  be  robust  in  the  presence  of  significant  noise.  The 
technique  is  intended  ultimately  be  combined  with  analytical  measurements  (e.g.,  Kalman  filter 
thermal/health  parameter  estimates)  for  FOD-event  detection  via  information  fusion  from  these 
(and  perhaps  other)  sources.  It  is  envisioned  that  fusion  of  structural  and  thermal  performance 
information  will  provide  a  pilot,  or  an  autonomous  vehicle,  more  confidence  in  a  FOD-event 
diagnosis  thus,  enabling  a  corrective  action  to  be  based  on  a  more  informed  decision. 
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Due  to  the  lack  of  high-frequency  FOD-event  test  data  in  the  open  literature,  a  reduced-order 
turbofan  structural  model  (ROM)  was  synthesized  from  a  finite  element  model  to  support  the 
investigation.  The  wavelet-based  FOD-event  detection  scheme  developed  is  “data  driven”  (i.e.,  not 
model-based),  thus  precise  correspondence  of  the  ROM  to  the  FEM  (which  itself  would  probably 
not  represent  an  actual  engine  with  complete  fidelity)  is  not  essential.  The  ROM  is  considered  to  be 
“phenomenologically  accurate,”  and  shown  to  provide  acceptable  fidelity  with  regards  to  the 
application  i.e.,  real-time  diagnostic  system  development. 

Using  the  state-space  ROM,  a  sensor  placement  analysis  was  performed  to  determine  the 
optimal  accelerometer  location  for  FOD-event  detection.  To  demonstrate  the  robustness  of  the 
technique,  testing  was  performed  using  the  least  optimal  bearing  location.  For  comparison,  analysis 
of  noise-free  and  noisy  accelerometer  signals  was  performed  using  a  standard  wavelet 
decomposition  and  shown  to  be  inadequate  for  detecting  FOD  events.  In  the  presence  of  significant 
noise  (SNR=3.5,  1 1  dB),  precise  location  of  the  FOD  event  in  time  was  obtained  using  the  custom 
wavelet  feature  developed.  Choice  of  the  smoothing  wavelet  used  (  Daubechies  8)  was  subjective, 
and  based  on  knowledge  of  the  signal  characteristic  of  interest  i.e.,  a  FOD  event-induced  abrupt 
change  in  an  otherwise  noisy  signal. 

Future  work  should  include  determining  the  optimal  wavelet  (and  corresponding  scale)  to  apply 
based  on  current  information  theoretic  techniques  such  as  maximum  entropy.  Additionally, 
combination  of  structural  and  thermal  performance  information  could  provide  a  robust  diagnosis-of 
foreign  object  damage  in  specific  engine  components  i.e.,  extending  this  research  to  the  Fault 
Detection  Isolation  (FDI)  problem.  Continued  development  and  verification  of  the  feature 
presented  would  require  high-frequency  accelerometer  data  (i.e.,  sensor  bandwidth  greater  than  10 
kHz)  from  FOD  event  engine  tests  at  multiple  operating  conditions. 

6.0  Nomenclature 

Am  Reconstructed  approximation  component  of  wavelet-analyzed  signal  at  scale  m 

M  Mass  Matrix 

C  Damping  Coefficient  Matrix 

Dm  Reconstructed  detail  component  of  wavelet-analyzed  signal  at  scale  m 
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DTWTD 


DTWTA 

K 

G 

Cxy,Kxy 

a,p 

co 

F,f 

F 

l0>  la 
m 

n 

r|(t) 

Y(f) 

Vm,n 

q(t) 


Discrete-time  wavelet  transform  detail  coefficient 

Discrete-time  wavelet  transform  approximation  coefficient 

Stiffness,  Stiffness  Matrix 
Gyroscopic  Matrix 

Quantities  with  coupling  between,  e.g.,  the  x  and  y  degrees  of  freedom 
Proportional  damping  coefficients 
angular  frequency  (rad/sec) 

Frequency  (Hz) 

Force 

Moments  of  inertia  about  the  X  and  Y  axes  respectively 
Dyadic  scale  of  wavelet-decomposed  signal 

Wavelet  time  delay  in  number  of  samples 

Displacement  in  modal  coordinates 
Discrete  Fourier  Transform  of  y(n) 

Wavelet  at  scale  m  and  translation  n 

Displacement  in  physical  coordinates 

Wavelet  scaling  function  at  scale  m  and  translation  n 
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